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Abstract 



Magnetohydrodynamic (MHD) and two-fluid quasi-neutral equilibria with azimuthal 
symmetry, gravity and arbitrary ratios of (nonrelativistic) flow speed to acoustic and 
Alfven speeds are investigated. In the two-fluid case, the mass ratio of the two species 
is arbitrary, and the analysis is therefore applicable to electron-positron plasmas. The 
methods of derivation can be extended in an obvious manner to several charged species. 
Generalized Grad-Shafranov equations, describing the equilibrium magnetic field, are 
derived. Flux function equations and Bernoulli relations for each species, together 
with Poisson's equation for the gravitational potential, complete the set of equations 
required to determine the equilibrium. These are straightforward to solve numerically. 
The two-fluid system, unlike the MHD system, is shown to be free of singularities. It 
is demonstrated analytically that there exists a class of incompressible MHD equilib- 
ria with magnetic field-aligned flow. A special sub-class first identified by S. Chan- 
drasekhar, in which the flow speed is everywhere equal to the local Alfven speed, is 
compatible with virtually any azimuthally symmetric magnetic configuration. Poten- 
tial applications of this analysis include extragalactic and stellar jets, accretion disks, 
and plasma structures associated with active late-type stars. 
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1 Introduction 



Observations across a range of wavelengths dating back over several decades have es- 
tablished the existence of jets associated with stars (Camenzind 1998; Henriksen 1998) 
and active galactic nuclei (Zensus 1997; Ferrari 1998). These structures, and also 
accretion disks, have in common the following features: they contain flows of bulk 
matter; they are embedded in magnetic fields; and they can be treated, to a first 
approximation, as having an axis of symmetry. All three features are also found in 
most of the laboratory magnetic confinement systems which have been investigated 
experimentally, such as tokamaks (Wesson 1997). There is increasing evidence gen- 
erally that magnetic fields play an important role in governing the physics of a wide 
variety of astrophysical processes. These include, for example, the formation and early 
evolution of galaxies, and interstellar gas dynamics (Zweibel & Heiles 1997). In the 
case of accretion disks, it is believed that purely hydrodynamic models cannot provide 
the effective viscosity required to account for observationally-inferred accretion rates, 
and that magnetohydrodynamic (MHD) effects must be taken into account (Hawley & 
Stone 1998). Magnetic fields may also provide a means of connecting accretion disks 
with jets, via the generation of MHD waves (Tagger & Pellat 1999). Finally, evidence 
has emerged recently that plasma confinement by dipole-like magnetic fields can ac- 
count for X-ray and radio emission from active late-type stars (Kellet, Bingham & 
Tsikoudi 2000). 

The first step in the construction of a theoretical model of a quasi-steady magnetized 
plasma structure (whether astrophysical or in the laboratory) is a determination of 
its equilibrium state. In general, this requires all time derivatives to be set equal to 
zero in the combined system of Maxwell and magnetized fluid equations, and solutions 
determined for the magnetic field, density, and, if applicable, flow velocity in three di- 
mensions. Having determined the equilibrium configuration, observed time variations 
can then often be interpreted theoretically as perturbations of the equilibrium state. 
MHD equilibrium studies of axisymmetric systems with flow have been carried out by 
many authors, including Chandrasekhar (1956), Woltjer (1959), Morozov and Solov'ev 
(1963), Zehrfeld & Green (1972), Maschke & Perrin (1980) and Throumoulopoulos k 
Pantis (1989). Early work in this field was restricted to the case of incompressible flow 
(Chandrasekhar 1956; Woltjer 1959). The problem of MHD equilibrium in toroidal 
systems with compressible flow was studied in the ideal limit by Morozov & Solov'ev 
(1963), and in the resistive case by Zehrfeld & Green (1972). The case of purely toroidal 
flow has been studied by Maschke & Perrin (1980) and Throumoulopoulos & Pantis 
(1989). Self-similar flow solutions were obtained by Blandford & Payne (1982), and 
Lovelace et al. (1986) developed a general theory of axisymmetric MHD equilibria 
with relativistic flows, which they applied to disks associated with rotating magnetized 
stars and black holes. Rosso & Pelletier (1994) used a variational method to resolve 
mathematical problems arising from MHD flow singularities. Bogoyavlenskij (2000) re- 
cently demonstrated the existence of exact axisymmetric MHD equilibria which do not 
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include flow, but may nevertheless have applications to astrophysical jets. Krashenin- 
nikov et al. (2000) included flow effects in a study of magnetic dipole equilibria, while 
Keppens & Goedbloed (2000) examined axisymmetric stellar wind equilibria with both 
open and closed magnetic field regions. 

These analyses of axisymmetric MHD equilibria are generally based on the Grad- 
Shafranov equation, derived originally to describe magnetic field equilibria in nuclear 
fusion experiments (Shafranov 1958; Grad & Rubin 1958): it is essentially an expression 
of momentum balance for one or more magnetized fluids. The analysis in this paper 
is also based on various forms of the Grad-Shafranov equation: we investigate ideal 
axisymmetric MHD (Section 2) and two-fluid (Section 3) equilibria with flows which 
are nonrelativistic but otherwise arbitrary. Whereas MHD equilibria have been studied 
in considerable detail by previous authors (in particular, Lovelace et al. 1986), little 
attention has been paid to two-fluid effects. We use a similar formalism for the MHD 
and two-fluid models, thus making it straightforward (and instructive) to compare 
and contrast them. The two-fluid model provides the basis for a more comprehensive 
description of jets and accretion disks than the MHD model, and is actually more 
tractable numerically. In the MHD case, special classes of solutions can be identified 
analytically (Section 4), which provide useful benchmarks for more realistic numerical 
solutions. 



2 General equilibrium analysis: MHD 

We present an alternative derivation of the "generalised Grad-Shafranov" equation of 
ideal MHD with arbitrary flows in azimuthally symmetric systems. This equation was 
first obtained by Lovelace et al. (1986): our alternative derivation is simpler than 
that of Lovelace et al., and can be readily generalized to the two-fluid case. A similar 
analysis was carried out by Goedbloed & Lifschitz (1997) for a system with translational 
rather than azimuthal symmetry. We consider non-relativistic MHD equilibria of a 
quasi-neutral plasma. The crucial simplifications are due to the assumed azimuthal 
symmetry (about the z-axis) and steady conditions. We adopt a cylindrical coordinate 
system, denoting the azimuthal angle by 0, and distance from the symmetry axis by 
r. It is useful to consider an "external" source of gravitation, creating an azimuthally 
symmetric gravitational potential V(r,z), although it will be seen that Newtonian 
self-gravitation of the plasma can be easily incorporated into the analysis. 

The Maxwell equation V • B = ensures that the magnetic field B has potential 
representation 

1.9* „ 1.9* 



B 



-e r + H — e 2 

r oz r or 



In equation (1) \l/(r, z) is the poloidal magnetic flux function and B^r, z) is the toroidal 
field. It is easily shown that \l//r is the toroidal (azimuthal) component of the magnetic 
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vector potential A. Since B • = 0, the function \l> is constant along magnetic 
field lines. In an analogous manner, mass continuity ensures that the plasma mass 
flux vector, defined to be the product of mass density p and flow velocity v, can be 
represented in the form 



pv 



Idx 1 dx 

r oz r or 



(2) 



where x is the mass flow function. Because of the assumed azimuthal symmetry, this 
is simply related to the poloidal magnetic flux function \l/ as follows. The ideal MHD 
Ohm's law in Gaussian cgs units is 

V<3> = — x B, (3) 

c 

where c is the speed of light and $ has gradient equal to minus the electric field E: 
the existence of such a potential follows from the steady-state assumption. Since $ 
cannot depend on 0, the azimuthal component of equation (3) yields 

v r B z = v z B r . 

Expressing the velocity and magnetic field components in terms of ^ and x using 
equations (1) and (2), we obtain 

d(x, *) _ Q 
d(r, z) 

from which it follows that 

X = F(n (4) 

where F is an arbitrary function. This, and other arbitrary functions appearing in 
the MHD and two-fluid systems of equations, are determined ultimately by plasma 
transport processes (Freidberg 1982). However, valuable physical insights can often be 
gained by adopting simple forms for these functions (see, e.g., Bogoyavlenskij 2000). 

The r and z components of equation (3) can be written in the form 

dr r dr ' 

and 

dz r dz' 

where F' = dF/d^. Introducing the quantity Q = {v$ — B^F' / p)/r, and eliminating 
$ in the two equations above by cross-differentiation and subtraction, we find that 

ggVg) 

<9(r, z) 
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This indicates that fHs a function of *. From the definition of fi, we thus have 



(7) 



In a similar fashion, by eliminating Q from equations (5) and (6), it is straightforward 
to show that $ is also a function of \1/ and that c$' = Q (the prime again denoting 
differentiation with respect to \I/). These results have several interesting physical in- 
terpretations. First, azimuthally symmetric, steady MHD equilibria with flows have 
electrostatic potentials which are functions of the poloidal magnetic flux, regardless of 
centrifugal and Coriolis effects at arbitrary flow Mach number. Second, in any poloidal 
plane (i.e. any (r, z) plane with <fi fixed), although not necessarily in three dimensional 
space, the flow and magnetic field components are parallel. Third, the poloidal mass 
flow function depends only on the poloidal magnetic flux. However, in general the 
density is not a flux function. 

To proceed further, it is necessary to calculate the components of the vorticity K = 
V x v, the current density j = (c/47r)V x B, and the vector products (V x B) x B, 
K x v. The following relations are easily derived from the definitions: 



V x B 



dB s 19,, 
-e r + j^e + - — {rB^e, 



dz 



r dr 



(8) 



J</> 



d 2 ^ d AdV s 
h r — ( 

dz 2 dr r dr ' 



K 



dv^ 



dz 



e r + K1e d 



ld_ 

r dr 



[rvd>)e. 



K l = — r 



d F' dtf +r d ^F'd$>. 

dz p dz dr rp dr ' 



(V x B) x B = 



. .,19* rBs d , n ,, 



r dr 



dr 



+(;; 



ld^ 



1 9(*,rS^) 
r 2 d(r, z) 

rBd d . 

\rB c 



r dz r 2 dz 



K x v 



F'dV r V<p d F> d^,r V4> ) 

( K <t>—^> 2 7r( rv <t>)) e r + —^7 f-e^ 

v pr dr r z dr pr z d{r,z) 

^ ^ pr dz r 2 dz 
The quantity j% — ^j^/c, where is the azimuthal current density. 



(9) 
(10) 

(11) 



:i2) 



(13) 



The system of MHD equations is completed by the ideal isentropic equation and the 
equation of motion. For a perfect gas with ratio of specific heats 7 the specific entropy 
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s is a function of p/ p 7 = a. In the absence of dissipation, a must be conserved in the 
fluid frame: 

v • V<7 = 0. 

It is clear from equations (2) and (4) that this condition is equivalent to 



d(r, z) 



0. (14) 



The steady-state equation of motion is the usual Eulerian one, with Lorentz, pressure 
and gravitational forces included. In terms of vorticity K it can be written in the form 

1 v 2 1 

Kxv = — Vp-V— -W+- — (VxB)xB. (15) 
p 2 4np 

To obtain equation (15) we have used the vector identity 

v 2 

v Vv = (V x v) x v + V— . (16) 

Equation (14) has the immediate consequence that a (a measure of the specific entropy) 
must be an arbitrarily specifiable flux function. 

We next consider the azimuthal component of the equation of motion. Noting the fact 
that all quantities are independent of 0, due to the assumed azimuthal symmetry, we 
infer from equations (12), (13) and (15) that 

F'd(^,r V(j) ) 1 <9(^,r^) 
r 2 d(r, z) Anr 2 d(r, z) 

Cancelling r 2 in the denominator on both sides and observing that F' depends only 
upon \1>, we see that this relation is equivalent to 



d(F»rv+ - IjB^ jf) 
<9(r, z) 



= 0. (17) 



This relation, an expression of the law of conservation of canonical angular momentum 
of the fluid about the symmetry axis, connects the toroidal flow and magnetic field 
components. Defining 

k = F'rv 4> -^-rB 4n (18) 
it is clear from equation (17) that A is a flux function. 

It will be seen that the quantities defined by equations (7) and (18) play a key role in 
the reduction of the equations of motion. Making use of equations (12) and (13), we 
find that equation (15) in the (r, z) plane can be written in the form 



w - in 



W 1 



V 



1 



pr 2r 2 L J 87rpr 2 

7 



1, 



rV ~t>) ~Q 2 V ( r ^) "-VP 



.v 2 



"Vy - W. (19) 

If the gravitational potential V and the "structure functions", F(ty), fi(\l/), A(\&) and 
cr^) are prescribed, equation (19) represents two partial differential equations for the 
determination of \I/ and p as functions of r, z, subject to suitable boundary conditions. 
The explicit forms of these equations will now be derived in full generality. 

We note that the left hand side of equation (19) is annihilated by taking the scalar 
product with the vector V\I/ x e^, which is always in a poloidal plane and tangential to 
the flux curves defined by constant values of ^f. The right hand side is simply related 
to the tangential derivatives of various quantities, by virtue of a geometrical relation 
valid for any function /(r, z): 

( v * xe *w = ^# = i v *if- < 20 > 

where / denotes arc length along the flux line. Using this, we annihilate the left hand 
side of equation (19) and derive the relation 



1 d[(rv*)*\ 1 d[(rB^] _ d \ 7 x + + y 



2r 2 dl 8vrpr 2 dl dl 



7 - r ' 2 



f21^ 



This equation represents a far-reaching generalization of the well-known Bernoulli 
equation of gas dynamics. A formal integral of it will now be obtained. 

Taking account of equation (18), we introduce the following representations of rv^ and 

nv = J(tt) + 6(r,z), (22) 
rB^ = /(*) + 47rF / (^)0(r, z). (23) 
Equation (18) can then be written in the form 

A(«) = JF' - £. (24) 

In equations (22) and (23) O is a function to be determined and /, J are arbitrary 
functions of \P related to A and Q. Substitution of these equations into equation (7) 
yields the relation 

(J + 0) - y(/ + 47rF , 0) =r 2 Q(^). (25) 

This equation expresses in terms of the flux functions Q, I and J, and in terms 
of p and r 2 , which are not necessarily flux functions. The quantity measures the 
variation of rv^ and rB^ on flux surfaces. In terms of J and 0, the two terms on the 
the left hand side of equation (21) are given by 

i Hr^_ i 9 I(j!+we + e!) , Ip + e)] ae 



2r 2 dl 2r 2 dl 



dl ' 
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8npr 2 



dl 



8irpr 2 dl 
1 

r 2 



(J 2 + 2/47rF'0 + 167r 2 (F') 2 e 2 ) 
dO 



IF' + 47r(F') 2 e 



dl 



Thus we have 
2r 2 dl 



d \{rB 4 



87ipr 2 



dl 



;j + 0) - y(/ + 47rF / 0) 



dl dl [ Uh 



(26) 

where we have made use of equation (25) and the fact that fl is a flux function (i.e. 
80. /dl = 0). Equation (21) can now be exactly integrated by writing it in the form 



d_ 
dl 



OQ 



7 



7 



■{op 



7-l> 



V 

T 



V 



0. 



From equation (25) it follows that 9 is given by 

r 2 + (yl - J) 







1 - 



(27) 



(28) 



It is clear from equation (27) that there exists a generalized Bernoulli integral 

ne '—(ap^- L ) - — -V = -h{V), (29) 



where the arbitrary function h{^) may be regarded as an MHD generalization of the 
"stagnation enthalpy" (total pressure) of gas dynamics (Meyer 1971). In the present 
calculation it is convenient to use the flux function 



= h(y)-n(y)j{v), 

rather than h itself. The Bernoulli relation then becomes 

2 

^ + V = if (tt) + Orv^ - —L- ap^-\ 
I 7 — 1 



(30) 



(31) 



Using the Bernoulli relation in the form given by equation (31) to eliminate v 2 /2 + V, 
we find that equation (19) can be written in the form 



k ; f ' - in 



(rB 4 



-A' + 



{rB^)(rv 4 



-F" + prH' + pr(rv<t)Q! 



rp a 



7 — 1 a 



W = 0. 

(32) 

Nontrivial solutions of this equation have V 1 !/ 7^ 0: setting the bracketed quantity in 
equation (32) equal to zero, and using our expressions for K£, j% [equations (9) and 
(11)], we infer that \I/ and p satisfy 



d 2 ^ d ,\d^, 
V r — ( 

dz 2 dr r dr ' 



AttF' 



d F'dV d F'dV. 
dz p dz dr rp dr 
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Atit 2 

+47rpr 2 H' + A-np^Vt'rv^ a' p 1 - AnrB^A' + ^(rB^rv^F" = 0. 

7-1 

This is a generalized form of the Grad-Shafranov equation. It is convenient to rearrange 
it in the form 

d dty d Id* C\I> 2 + $ 2 ) 

+ + 4vrF'F" ^ + WJ 

oz oz or r or p 

4:7rr^ 

+4npr 2 H' + A-npr^Vt'rv^ & ' p< - AnrB^A' + Air^B^) (rv^F" = 0, (33) 

7-1 

where 

A(p,F') = l-^p^. (34) 

It is straightforward to verify that Eq. (33) is identical to the Grad-Shafranov equation 
obtained by Lovelace et al. (1986). Given the arbitrary flux functions F, A, Q, H, 
a and the potential V, equations (7), (18), (31), (33) and (34) determine * and p. 
The coefficients of the highest order derivatives in equation (33) vanish, and hence the 
equation becomes singular, when A = 0: physically, this corresponds to the projection 
of the flow velocity onto any (r, z) surface being equal to an Alfven speed defined 
in terms of the (r, z) components of B. The fact that the density p depends, via 
the Bernoulli relation [equation (31)], on * and its derivatives means that equation 
(33) contains other singularities (Lovelace et al. 1986), which are less immediately 
apparent than the one corresponding to A = 0. The presence of p in equation (33) 
also means that it is fundamentally nonlinear, regardless of the choice of arbitrary 
functions: in this respect it differs from the Grad-Shafranov equation without flows 
(e.g. Bogoyavlenskiy 2000). The existence of flow singularities aggravates considerably 
the difficulties involved in finding solutions (see e.g. Rosso & Pelletier 1994). In 
principle, the singularities can be removed by invoking dissipation or electron inertia, 
or by special, compatible choices of the arbitrary functions involved. 

We have so far assumed that the gravitational potential V is due to a distribution of 
masses which is external to the plasma. It is straightforward to include Newtonian 
self-gravitation of the plasma by adding the gravitational Poisson equation 

Id dV d 2 V 

;*( r *' + a?= 4 ^" + "-'- < 35) 

where p ext represents the density distribution of external, uncharged bodies (e.g. a 
neutron star or black hole). The Poisson integral solution of this equation can be used 
to eliminate V from the generalized Bernoulli relation [equation (31)]: this would have 
the effect of making the governing equations non-local. 
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3 General equilibrium analysis: two-fluid theory 



It is possible to generalize the above model to include two-fluid effects. There are both 
physical and mathematical reasons for seeking such a generalization. Ferrari (1998) has 
noted that single-fluid models are unlikely to describe adequately the microphysics of 
extragalactic jets, and that a two-fluid theory would provide a useful intermediate step 
towards the development of a fully kinetic model. Moreover, the two-fluid equations of 
motion include inertial terms which, as noted above, remove troublesome singularities 
appearing in the MHD theory. It is convenient to have a symmetrical notation for ions 
(mass mj, charge = e) and electrons (mass m e , charge e e = — e). In principle, the 
ions can be any charged species, including positrons. Indeed, the methods are also 
applicable, mutatis mutandis to quasi-neutral plasmas with several charged species. 
The equations derived below are non-relativistic, assume quasi-neutrality, azimuthal 
symmetry and steady conditions. Electron inertia and temperature will be included, 
but all dissipative and irreversible terms are neglected. 

As before, we use the magnetic field representation given by equation (1). In place of 
the mass density p we introduce the common number density, n = n e = rij. We replace 
equation (2) by equations for the particle flux functions, the existence of which follows 
from the continuity equations for the two species (xj ~ * @j m j> 3 = h e ) 



1 90, 1 <99 7 

• — 7T- e r + nvtfej H ^— e 2 

r oz r or 



(36) 



where Qj, v^j are, respectively, the particle flux function and toroidal flow speed of 
species j. Introducing the vorticities Kj of the two species, we infer a set of relations 
analoguous to equations (10), (11) and (13): 



K, 



dz 



19, , 



(37) 



<t>3 r Q r 



dz n dz 



d , 1 



dr rn dr 



rnv 



4>:i 



d 



dr 



1 d(<d h r V<t)j ) 



d(r, z) 



e<f, 



IdQj rnv ( p j d 



4>3 r dz 



dz 



{rv^))e z . 



(38) 



(39) 



The two equations of motion governing momentum balance now take the form 



Kj x nVj = —Vpj - mjnV\ 2 J2 - m^VV - e^nV^ + ejnVj x B/c, (40) 



where pj is the pressure of species j, = e = — e e and, as before, $, V are, respectively, 
electrostatic potential and gravitational potential. Adding the electron and ion equa- 
tions of motion [equation (40)], we obtain the single-fluid MHD equation of motion 
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[equation (15)] in the formal (singular) limit e — > oo (keeping all other quantities fixed). 
This corresponds to the ion Larmor radius and collisionless skin depth [c divided by 
the electron plasma frequency uj pe ] both tending to zero. In this limit, the electron 
equation of motion reduces to the ideal MHD Ohm's law [equation (3)]. 

We begin by considering the azimuthal component of equation (40). Substituting from 
equation (39) and evaluating the azimuthal component of Vj x B, we obtain 



d(r, z) 



0, 



where IL,- are canonical momenta of the two fluids [cf. equation (18)]: 

rif = mjrv^ + ej^/c, 



(39) 



(42) 



whence it follows that 



(43) 



where Fj are arbitrary functions of the respective particle flux functions. We now apply 
Ampere's law, 

Ane 

V x B = (nvj - nv e ). (44) 

c 

Only two of the three components of this are independent. The r and z components 
integrate to give 



rB 4> = — (Qi-Qe), 
while equation (8) indicates that the component can be written as 

Alien 



(45) 



(46) 



Using equations (9), (42) and (43) we rewrite equation (46) in the form 



d 2 ^ 9,19*. 
1_ r — ( ) 

dz 2 dr r dr 



Anen 



El 

rrii 



eirrii + m e )^ 



This is the Grad-Shafranov equation for the two-fluid system. 



(47) 



We assume, in the absence of dissipation, that the entropies of the two species are 
constant following the flow: 



dl 



0, 



(48) 



where Oj = Pj/ri 1 . Taking dot products of Vj with equation (40), and using equation 
(47), we derive two Bernoulli relations 



7 — 1 n Z 



(49) 
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where, as in the MHD case, the stagnation enthalpies Hj(Qj) are arbitrarily prescrib- 
able functions of the respective particle flux functions. We note that these two equa- 
tions effectively determine the number density n and the electrostatic potential $, 
given <7j(0j), Fj(Qj), Qj, * and V. Since <3> appears nowhere else explicitly, it may be 
eliminated by simply adding the two equations to obtain the following equation for n: 



+ a e )n^ + + ^) + K + m e )V = H t (Q t ) + H e (Q e ). 

7 — 1 Z Z 

Using equation (49) we may rewrite equation (40) in the form 



(50) 



e,- 



m. 



Kj + -^B) x (nvj) = -nH'VQj + 



a 



c J 7 — 1 (Xj 

From the radial or z components of this we immediately obtain 



(51) 



F'AFj - ^) = -nrH' + 



7 - 1 <Tj ' 



i.e. 
3 199 



9 , 1 96, 



9z n 9z 



9r rn dr 



ej „ n „, , „ e 



rB*+-*F;(F j —±*) = nr 2 H' 3 /m r 



2r»7 



r n 



m 



771,(7- 1) 
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Eliminating r_B^ using equation (45), we obtain the following closed system for Qj, *, 
n and V: 



9 ,1 90i, 9 , 1 96^ 

— ( -) + r — ( ~) 

dz n dz dr rn dr 



' lTC .( e< - e e ) + -VOT 



TTljC C 



= nr 



..2 i 



r 2 n 7 



mi 77^(7 - 1) 



mt 



9 ,190 e , 9 , 1 99 e / 

— ( -) + r — ( -) 

dz n dz dr rn dr 



+ — — (Qi - 8e) + + 

m e c c c 



= nr 



r 2 n 7 



m e m e (7 - 1) 



7 



9 2 * 9,19*/ 
h r — ( ) 

9z 2 9r r 9r 

,7-1 
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m,-v 2 m R v 



-fo + a e )tf- 1 + (^ + ^p) + k + m e )y = ^(eo + # e (e e ), 

1 9 , 9F S 9 2 1/ 
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-(r— ) + 

r 9r 9r 9z 2 



47rG [(m^ + m e )n + p cxt ] . 



(52) 



(53) 
(54) 
(55) 
(56) 
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The first two equations are the equations of motion, the third represents Ampere's law, 
equation (55) is the Bernoulli relation for the total pressure, and equation (56) is the 
gravitational Poisson equation. 

Given the six arbitrary functions Fj, <jj, Hj and the "external" mass density p ext , these 
equations have to be solved subject to suitable boundary conditions on 0j and 
the gravitational potential V. They fully determine the structure of both flows and 
magnetic fields. The velocities, currents and so on can be obtained from the various 
auxiliary relations. This system represents an exact generalization of the equations 
derived by Lovelace et al. (1986) to two-fluid, dissipationless, azimuthally symmetric, 
quasi-neutral, nonrelativistic, gravitating equilibria with arbitrary flows. It is impor- 
tant to note that the leading order operators in the partial differential equations are 
entirely nonsingular and elliptic, provided the number density (n) remains bounded. 
For this reason, the system is much easier to deal with numerically than the equations 
of MHD equilibrium. 



4 Special solutions: field-aligned flows 

In this section a particular class of MHD equilibria is considered. In Section 2 it was 
shown quite generally that the flow function x is a function of the magnetic flux, \I>: we 
now assume that mass density p is also a flux function. In this case, it is easily shown 
from equations (2) and (4) that the flow velocity is divergence-free, i.e. the plasma is 
incompressible: in general, this requires that the flow speed v be less than the sound 
speed c s (e.g. Landau & Lifshitz 1987). Axisymmetric MHD equilibrium equations 
were obtained for general incompressible flow by Woltjer (1959). Here we consider 
incompressible flows with fl = [cf. equation (7)]: in such cases rv^ is proportional to 
rB^. Denoting F'/p by G' [cf. equations (2) and (4)], we obtain 

v = G'B, (57) 

Where G($f) is a new flux function. This equation indicates that the flow and the field 
are aligned in three dimensions [in general, the alignment is in (r, z) planes only]. The 
ideal Ohm's law [equation (3)] indicates that $ = for such flows. 

From equation (57) we infer an expression for the vorticity K = V x v in terms of G 
and B: 

K = G"'W xB + G'VxB. 
Substituting this in the equation of motion [equation (15)], we get 

v 2 1 

pG' [G"'W x B + G'V x B] x B = -Vp - pV pW + — (V x B) x B. 

2 47T 

Using B • V 1 !/ = 0, we find that the above equation may be written in the form 

1 v 2 
-pG"G"'B 2 W + (pG" 2 - — )(V x B) x B = -Vp - pV— - pW. (58) 

4"7T 2 
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If we take the dot product of this equation with B, consistency requires that the right 
hand side must vanish. Using the fact that p is a flux function, we infer from this that 



v 2 



B • V(p + p— + pV) = 0. 



We thus obtain the Bernoulli relation 

P + PY + PV = H*W, (59) 

where H* is an arbitrary function of Using this to eliminate p from equation (58), 
we obtain 

-pG"G"'B 2 W + (pG' 2 -^)(VxB)xB = -iT'W + p'W(y + V). 
Putting v = G"B [equation (57)], this reduces to 

-{pG'G" + p'^)B 2 W + ( P G' 2 -^)(VxB)xB = -(#*' - pV)W. (60) 

This is our fundamental equation. We first note that since \1/ depends only on r and 
z, the azimuthal component is satisfied if j x B vanishes in the <fi direction. It is clear 
from equation (12) that this requires rB^ to be a flux function, which we identify with 
7(\1/) defined by equation (23) [G = in this case]. It follows from our field alignment 
assumption that rv$ = J(^) = G'l. We observe next that if pG' 2 is a constant, the 
first term on the left hand side of equation (58) vanishes. In this case, making use of 
the r or z components of equation (12), we find that the equation reduces to 

m 1 - - ^) = -(**'- p'n (6i) 

where 

Equation (61) can be simplified, using equation (9), to give a Grad-Shafranov equation 
for field-aligned MHD flows: 



A 



d 2 ^ d AdV, TTI 
1- r — ( ) + //' 

dz 2 dr r dr 



+ r 2 p'V = r 2 H*'. (62) 



This is algebraically similar to the Grad-Shafranov equation used by Bogoyavlenskiy 
(2000), but represents a generalization of it to include field-aligned flow in a spatially- 
varying gravitational potential. From the definition of G' [equation (57)], it is clear 
that 4npG' 2 = v 2 /c 2 A , where ca = B / 'y/Airp is the Alfven speed. Although we assumed 
that pG' 2 was a constant in order to obtain equation (62), this does not, of course, 
preclude the possibility of v, B and p individually varying in space. 
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There are several limiting cases of equation (62) which are of physical interest. If V is 
a prescribed gravitational potential, p and H* are quadratic in and IV is linear in 
we obtain a linear eigenvalue problem for \1/ and A, which can be solved by a variety 
of methods (e.g. Bogoyavlenskij 2000). Even when equation (60) is nonlinear, it can 
be solved numerically in a straightforward way. 

A special case of great interest is that of constant density and force-free (Beltrami) 
magnetic fields, defined by the property V x B = /iB for some scalar function /i, so 
that (V x B) x B = (Lundquist 1951; Mahajan & Yoshida 1998). Consistency with 
equation (60), with constant pG' 2 ^ 1/47T, requires that H* also be a constant. More 
generally, the force-free condition is satisfied for constant pG' 2 whenever H*' = p'V. 
From equation (62), it is clear that A is then an arbitrary finite constant, and that \l/ 
satisfies the equation 

d 2 ^ d 1 d^f 

Having solved equation (63) for one can obtain the velocity field v(r, z) from equa- 
tion (57): the pressure is then determined by the Bernoulli relation [equation (59)], as 
befits an incompressible flow. 

A second class of solutions has A = 0. The equilibrium is then compatible with a 
completely arbitrary azimuthally symmetric field! It is clear from equation (60) that 
it is not necessary for the azimuthal component of j x B to vanish in this case: the 
only consistency requirements are that pG' 2 is a constant and that H*' = p'V . Since 
A = pG' 2 — 1/An = {v 2 /c\ — l)/47r, the condition A = means that the flow is 
everywhere "trans-Alfvenic" (v = ca), as well as being field-aligned. This solution 
was first obtained by Chandrasekhar (1956), who proved moreover that it is stable. 
Consistency with the assumption of incompressibility (V • v = 0) requires in this case 
that c s > ca, i.e. that the plasma beta (3 ~ c 2 /c\ be greater than unity. It is not clear 
how restrictive this condition is, since (3 in most astrophysical plasmas is not accurately 
known. 

The most general field-aligned flow governed by equation (59) is one in which A, and 
hence pG' 2 , are flux functions, the equation then taking the form 



where 



d 2 ^ d ,ldV. TT , 
1_ r — ( ) + //' 

dz 2 dr r dr 



+ iA'((W) 2 + / 2 ) =r 2 (H*' - p'V), (64) 



m = PG' 2 - ^- (65) 

It is plain that we can extend the same ideas to the two-fluid case by taking n to be 
a flux function and Vj = (^(^B. Specific applications of these solutions are left to 
future work. 
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5 Conclusions and Discussion 



We have derived equilibrium equations for azimuthally symmetric MHD and two-fluid 
systems with arbitrary nonrelativistic flows. Our method of derivation in the MHD 
case is distinct from that used by Lovelace et al. (1986), and provides the basis for 
our two-fluid calculation. We have identified some special exact solutions of the MHD 
system of equations which are of independent interest. Since there are no restrictions on 
the geometry (other than azimuthal symmetry), both systems of equations are equally 
applicable to astrophysical and laboratory systems with azimuthal symmetry. The 
two-fluid model is more suitable than the MHD model for astrophysical applications, 
not only because it contains more physics, but also from the computational point of 
view: the two-fluid equations do not involve singularities which are associated with 
ideal MHD (Blandford & Payne 1982; Lovelace et al. 1986). The physical reason for 
this is that the two-fluid model, unlike the MHD model, takes into account electron 
inertia, which introduces a new fundamental length, the collisionless skin depth c/uj, pe 
(in the MHD limit, c/uj pe — > 0). In future work we intend to solve numerically both the 
MHD and two-fluid systems of equations for specific astrophysical scenarios, including 
extragalactic jets. 

The equations we have derived are nonrelativistic. Flow speeds much less than c occur, 
for example, in jets, hot spots and lobes associated with radio galaxies, but close to the 
core v ~ c (Ferrari 1998). If the assumptions of stationarity and azimuthal symmetry 
are retained, it is reasonable to expect that similar but more complicated equations will 
exist for relativistic conditions (this was demonstrated by Lovelace et al. for the case of 
ideal MHD). Although we considered plasmas with only one or two separate species, the 
methods clearly extend to any number of charged species, provided that overall quasi- 
neutrality and other equilibrium assumptions apply. It is relatively straightforward to 
take into account higher order collisional or collisionless effects (using, for example, 
the reciprocals of the Reynolds or Lundquist numbers as expansion parameters): one 
would expect the inclusion of such effects to provide transport equations for the various 
arbitrary functions appearing in the MHD and two-fluid systems of equations. 
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